##__docformat__ = "restructuredtext en"
# ******NOTICE***************
# two_game.py module by Shoibal Chakravarty
#
# You may copy and use this module as you see fit with no
# guarantee implied provided you keep this notice in all copies.
# *****END NOTICE************
"""
This module contains the functions that simulate a scenario where two regions
have a two period planning horizon (current period and the next period, where each region is myopic), and the 
third region is myopic. The two regions play a Nash game and the solution is the intersection 
of the best response lines which is determined numerically. The green energy deployment in the 
last period is myopic.
"""

__author__ = "Shoibal Chakravarty (shoibalc@princeton.edu)"
__version__ = "0.1"
__date__ = "2011/August"
__license__ = "BSD"
__all__ = ['twomyopictwonash','twonashupdate','intersectionregion','findintersectionregion','finalperiodupdate']
import numpy as np
from utility import U2
from myopic import myopicupdate, regionupdate, worldupdate, cumulativecost
from flat2 import flatten
import scipy.optimize

def twomyopictwonash(r1,r2,r3,r4,w,com,precision=0.5,init_grid_size=40):
    """
    Function onemyopictwonash solves the scenario where regions r1 and r2
    deploy green energy by playing a two period Nash game (2nd period myopic)
    whereas r3 plays a myopic strategy. For every period except the last, r3 first
    plays the myopic strategy, and r1 and r2 play their best response Nash strategy.
    The best response is the green deployment that maximizes the region r1's two period utility for a
    given value of region r2's green energy deployment. This gives a best response curve as a function
    of r2's green energy deployment. This curve is calculated numerically for each region (r1 and r2).
    The solution to the Nash game is the intersection of the two curves which is also obtained numerically to 
    a given tolerance. In case of multiple solutions, a solution at random is picked.\n\n
    onemyopictwonash(r1,r2,r3,...):\n
    r1 and r2 play a two region Nash game with one period look ahead, r3 is myopic.
    """
    r3.descriptor = r4.descriptor = 'Myopic'
    r1.descriptor = r2.descriptor = 'Solution of a two region Nash game.'
    print "\n\nIn this scenario %s and %s  play a two region two period Nash game, %s and %s are myopic\n"% (r1.name, r2.name, r3.name, r4.name)
    for t in range(w.T-1):
#  handle the myopic region c first
        r3.g[t], r3.u1[t], r3.c1[t], r3.m[t] = myopicupdate(r3,r1,r2,r4,w,t)
        regionupdate(r3,t)
        r4.g[t], r4.u1[t], r4.c1[t], r4.m[t] = myopicupdate(r4,r1,r2,r3,w,t)
        regionupdate(r4,t)
# two region two period nash
        twonashupdate(r1,r2,r3,r4,w,com,t,precision,init_grid_size)
        worldupdate(r1,r2,r3,r4,w,com,t)
        print "\nPeriod", t+1,"\n", r1.name, "g", r1.g[t], "e", r1.e[t],"\n", r2.name, "g", r2.g[t], "e", r2.e[t],\
        "\n",r3.name, "g", r3.g[t], "e", r3.e[t], "\n",r4.name, "g", r4.g[t], "e", r4.e[t], "\ncg", w.cg[t+1],"\n"
# the last period is still myopic
    t = w.T-1
    finalperiodupdate(r1,t)
    finalperiodupdate(r2,t)
    r3.g[t], r3.u1[t], r3.c1[t], r3.m[t] = myopicupdate(r3,r1,r2,r4,w,t)
    regionupdate(r3,t)
    r4.g[t], r4.u1[t], r4.c1[t], r4.m[t] = myopicupdate(r4,r1,r2,r3,w,t)
    regionupdate(r4,t)
    worldupdate(r1,r2,r3,r4,w,com,t)
    print "\nPeriod", t+1,"\n", r1.name, "g", r1.g[t], "e", r1.e[t],"\n", r2.name, "g", r2.g[t], "e", r2.e[t],\
    "\n",r3.name, "g", r3.g[t], "e", r3.e[t], "\n",r4.name, "g", r4.g[t], "e", r4.e[t], "\ncg", w.cg[t+1],"\n"
    cumulativecost(r1,r2,r3,r4,w)


def finalperiodupdate(r,t):
    """
    Function finalperiodupdate writes the last period solution for the one period
    look ahead or Nash case. The last period solution is myopic and was calculated
    in the solution of the penultimate period when the two period problem was solved. This function
    copies the second period myopic solution of the penultimate period into the last period solution.
    """
    r.g[t] = r.g2[t-1]
    r.u1[t] = r.u2[t-1]
    r.c1[t] = r.c2[t-1]
    r.m[t] = min(r.M[t]*r.E[t+1],(r.E[t+1] - r.G[t]))
    regionupdate(r,t)


def twonashupdate(ra,rb,rc,rd,w,com,t,precision,init_grid_size):
    """
    Function twonashupdate solves the two region two period Nash problem for
    the current period (except for the last period in the simulation). The function
    successively calls function intersectionregion which return the current best guess 
    and a rectangle bounding the solution. This procedures continues till the desired 
    precision is reached.
    """
    maxa = min(ra.M[t]*ra.E[t+1],(ra.E[t+1] - ra.G[t]))
    mina = 1e-6
    maxb = min(rb.M[t]*rb.E[t+1],(rb.E[t+1] - rb.G[t]))
    minb = 1e-6
    D0 = 22
    D = D0
    EPS = EPSa = EPSb = init_grid_size*D
    print 'Search region',ra.name,'min = 0 max =', maxa, rb.name, 'min = 0, max =', maxb
    flag = True
    num = False
    xbest, ybest = np.nan, np.nan
    while EPS > precision and flag:
# initial run with  init_grid_size. Subsequent runs reduce it by a factor of D
        if num:
            D = np.random.randint(D0+1,D0+5)
        EPSa = EPSa/D
        EPSb = EPSb/D
        num = True
        print "GRID SIZE", D,
        mina, maxa, minb, maxb, xbest1, ybest1, lenxy = \
                intersectionregion(ra,rb,rc,rd,w,com,t,mina,maxa,minb,maxb,EPSa,EPSb)
        err_best = (xbest - xbest1)**2 + (ybest - ybest1)**2
        if np.sqrt(err_best) <= 0.25*precision:
# this hack prevents infinite looping if we hit a multiple equilibria. we select a particular solution and carry on.
# it also helps in avoiding an extra round of brute force maximization in some cases.
            if lenxy == 1:
                print "1. ATTAINED REQUIRED PRECISION"
                flag = False
            else:
                print "2. MULTIPLE SOLUTIONS DETECTED, PICKING AVERAGE OF DEGENERATE SOLUTIONS"
                flag = False
            break
        else:
            xbest, ybest = xbest1, ybest1
        EPSa = maxa - mina
        EPSb = maxb - minb
        EPS = max(EPSa, EPSb)
    if flag:
        print "3. ATTAINED REQUIRED PRECISION"
        ga = 0.5*(mina + maxa)
        gb = 0.5*(minb + maxb)
    else:
        ga = xbest
        gb = ybest
    
    parama = \
    (gb, com.cg0, com.greenpar, ra.B[t], ra.B[t+1], ra.cn[t], ra.co[t], ra.cn[t+1], ra.co[t+1], ra.G[t], rb.G[t], \
    rc.G[t]+rd.G[t], rc.G[t+1]+rd.G[t+1], w.G[0], ra.e[t], ra.E[t+1], ra.e[t+1], ra.E[t+2], ra.M[t], ra.M[t+1], (1-com.discount)**com.years,True )
    u, ra.g[t], ra.u1[t], ra.c1[t], ra.g2[t], ra.u2[t], ra.c2[t], ra.m[t] =  U2(ga,*parama)
    regionupdate(ra,t)
    
    paramb = \
    (ga, com.cg0, com.greenpar, rb.B[t], rb.B[t+1], rb.cn[t], rb.co[t], rb.cn[t+1], rb.co[t+1], rb.G[t], ra.G[t], \
    rc.G[t]+rd.G[t], rc.G[t+1]+rd.G[t+1], w.G[0], rb.e[t], rb.E[t+1], rb.e[t+1], rb.E[t+2], rb.M[t], rb.M[t+1], (1-com.discount)**com.years,True )
    u, rb.g[t], rb.u1[t], rb.c1[t], rb.g2[t], rb.u2[t], rb.c2[t], rb.m[t] =  U2(gb,*paramb)
    regionupdate(rb,t)


def intersectionregion(r1,r2,r3,r4,w,com,t,mina,maxa,minb,maxb,EPSa,EPSb):
    """
    Function intersectionregion takes among its parameters the coordinates of a rectangular  search region 
    and a grid-size for discretization.. It uses a brute force method (scipy.optimize.brute) to numerically 
    obtain the best response curves for regions r1 and r2. The best response curves are passed on
    to function findintersectionregion which returns the  intersection point(s) and the coordinates 
    of  smaller rectangular search region.
    """
#   Ga0,Gb0,Gc0,Gc1,G0,ea0,ea1,Ea0,Ea1,eb0,eb1,Eb0,Eb1,delta,precision,alpha, mina, maxa, minb, maxb, EPSa, EPSb = args
    lena = np.ceil((maxa - mina)/EPSa)
    lenb = np.ceil((maxb - minb)/EPSb)
    EPSa = (maxa - mina)/lena
    EPSb = (maxb - minb)/lenb
    eps = 1e-5
    slice_a = np.s_[mina: maxa + eps:EPSa,]
    slice_b = np.s_[minb: maxb + eps: EPSb,]
    mgrid_a = np.mgrid[slice_a[0]]
    mgrid_b = np.mgrid[slice_b[0]]
    listab = np.zeros((lenb+1))
    listaa = np.zeros((lenb+1))
    listbb = np.zeros((lena+1))
    listba = np.zeros((lena+1))
    m=0
    param =\
        [0,com.cg0,com.greenpar, r1.B[t], r1.B[t+1], r1.cn[t], r1.co[t], r1.cn[t+1], r1.co[t+1], r1.G[t], r2.G[t], \
        r3.G[t]+r4.G[t], r3.G[t+1]+r4.G[t+1], w.G[0], r1.e[t], r1.E[t+1], r1.e[t+1], r1.E[t+2], r1.M[t], r1.M[t+1], (1-com.discount)**com.years, False]
    for tt in mgrid_b:
        param[0]= tt
        parama = tuple(param)
        x0 = scipy.optimize.brute(U2,slice_a,args=parama,full_output=False,finish=None)
        listab[m] = tt
        listaa[m] = x0
        m+=1
#        parama = (tt, Ba, cfna, cfoa, Ga0, Gb0, G0, ea0, Ea0, ea1, Ea1, delta, alpha, True)
#        print tt, x0, U2(x0,*parama)
    m=0
    param = \
        [0,com.cg0,com.greenpar, r2.B[t], r2.B[t+1], r2.cn[t], r2.co[t], r2.cn[t+1], r2.co[t+1], r2.G[t], r1.G[t], \
        r3.G[t]+r4.G[t], r3.G[t+1]+r4.G[t+1], w.G[0], r2.e[t], r2.E[t+1], r2.e[t+1], r2.E[t+2], r2.M[t], r2.M[t+1], (1-com.discount)**com.years, False]
    for tt in mgrid_a:
        param[0] = tt
        paramb = tuple(param)
        #paramb = (tt, Bb, cg0, cfnb, cfob, Gb0, Ga0, Gc0, Gc1, G0, eb0, Eb0, eb1, Eb1, delta, alpha, False)
        x0 = scipy.optimize.brute(U2,slice_b,args=paramb,full_output=False,finish=None)
        listbb[m] = x0
        listba[m] = tt
        m+=1
#        paramb = (tt, Bb, cfnb, cfob, Gb0, Ga0, G0, eb0, Eb0, eb1, Eb1, delta, alpha, True)
#        print tt, x0, U2(x0,*paramb)
#    return listax, listay, listbx, listby
#    a = line(zip(listax,listay),color = 'blue')
#    a += line(zip(listbx,listby),color = 'red')
#    a.show()
    mina, maxa, minb, maxb, xbest, ybest, lenxy0   = findintersectionregion(listab,listaa, listbb, listba)
    print "Search region update:", r1.name, "min", mina, "max", maxa, "best", xbest, r2.name, "min", minb, "max", maxb, "best", ybest
#    if flag:
    return mina, maxa, minb, maxb, xbest, ybest, lenxy0


def findintersectionregion(x1,y1,x2,y2, Wa = 3, Wb = 3):
    """
    findintersectionregion finds the intersection point(s) of the two curves and returns 
    the intersection point and the coordinates of a  smaller bounding rectangular region. When there are more than
    one intersection points (i.e. degenerate solutions in this case) the algorithm return the coordinates
    of the rectangle bounding the intersection points and the  center of mass of these points. The center of mass
    of these points is also a solution as the degenerate solution space is this case is  a straight line.\n\n
    (Inspired by   intersections.m, function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust) of Douglas M. Schwartz.
    http://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections.)
    """
    eps = 1e-7
    n1 = len(x1) 
    n2 = len(x2)
    x1 = np.c_[x1]
    y1 = np.c_[y1]
    x2 = np.r_[x2]
    y2 = np.r_[y2]
# Find points that are in both lines.
    
    xb,xa = \
    np.nonzero(((np.tile(x1,(1,n2))- np.tile(x2,(n1,1)))==0.0)& \
    ((np.tile(y1,(1,n2))- np.tile(y2,(n1,1)))==0.0))
#    print xb,xa
    x0 = [x1[xx].item() for xx in xb]
    y0 = [y1[xx].item() for xx in xb]
    xy0 = zip(y0,x0)           # detects duplicates
    
# xbest0 and ybest 0 are useful if we have multiple solutions 
# (note that (xbest0,ybest0) is a solution in the degenerate set)
# instead of going into an infinite loop nextG2 checks if xbest0 
# and ybest0 change sufficeintly and breaks.
    xbest0 = np.mean(y0)
    ybest0 = np.mean(x0)
    print "CURRENT GUESS(ES)",
    print xy0
    xa1 = [xai - Wa > 0 and x2[xai - Wa].tolist() or x2[0].tolist() for xai in xa]
    ya1 = [xai - Wa > 0 and y2[xai - Wa].tolist() or y2[0].tolist() for xai in xa]
    xb1 = [xbi - Wb > 0 and x1[xbi - Wb].tolist() or x1[0].tolist() for xbi in xb]
    yb1 = [xbi - Wb > 0 and y1[xbi - Wb].tolist() or y1[0].tolist() for xbi in xb]
    xa2 = [xai + Wa < (n2 - 1) and x2[xai + Wa].tolist() or x2[-1].tolist() for xai in xa]
    ya2 = [xai + Wa < (n2 - 1) and y2[xai + Wa].tolist() or y2[-1].tolist() for xai in xa]
    xb2 = [xbi + Wb < (n1 - 1) and x1[xbi + Wb].tolist() or x1[-1].tolist() for xbi in xb]
    yb2 = [xbi + Wb < (n1 - 1) and y1[xbi + Wb].tolist() or y1[-1].tolist() for xbi in xb]
    mina, maxa, minb, maxb = min(flatten([ya1,ya2,yb1,yb2])), max(flatten([ya1,ya2,yb1,yb2])),\
        min(flatten([xa1,xb1,xa2,xb2])), max(flatten([xa1,xa2,xb1,xb2]))
    return mina, maxa, minb, maxb, xbest0, ybest0, len(xy0)


